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TRACTION BOUNDARY CONDITIONS FOR MOLECULAR STATIC SIMULATIONS 


XIANTAO LI AND JIANFENG LU 


Abstract. This paper presents a consistent approach to prescribe traction bound¬ 
ary conditions in atomistic models. Due to the typical multiple-neighbor interactions, 
finding an appropriate boundary condition that models a desired traction is a non¬ 
trivial task. We first present a one-dimensional example, which demonstrates how such 
boundary conditions can be formulated. We further analyze the stability, and derive its 
continuum limit. We also show how the boundary conditions can be extended to higher 
dimensions with an application to a dislocation dipole problem under shear stress. 


1. Introduction 

Atomistic models have established a critical role in material modeling and simula¬ 
tions. In order to study the mechanical responses, boundary conditions (BC) must be 
imposed. While specifying the displacement of the atoms at the boundary is straight¬ 
forward, imposing a traction is much more challenging due to the fact that the range of 
the atomic interactions typically goes beyond nearest neighbors. In direct contrast to 
continuum mechanics, where the boundary is of lower dimension (curves or surfaces), 
the ‘boundary’ in an atomistic model often consists of a few layers of atoms. As a result, 
there are multiple ways to prescribe a traction BC. For instance, forces can be applied 
directly to atoms at the boundary in such a way that they add up to the given traction. 
However, it is unclear how to distribute these forces among the atoms. In particular, 
boundary layers may develop and create large modeling error. 

Meanwhile, many mathematical problems associated with material defects have been 
formulated as a system under traction. Examples include cracks under mode-I loading 
(T^ , where uniform stress can be specified in the far field, and dislocations under shear 
stress Ij^, which led to the important concepts of Peierls barrier. Problems of this type 
can not simply be treated with BCs that prescribe the displacement of the atoms at the 
boundary Another possible approach to introduce traction is the Parrinello-Rahman 
method |[^ , where the stress is created by allowing the shape of the simulation cell to 
change, which is particularly useful when phase transformation processes occur. But 
the method is limited to periodic cells, and it can not treat material defects without 
introducing artificial images. 

The purpose of this paper is to formulate a proper BC that represents a traction force 
along the boundary. We set up the problem by embedding the computational domain 
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within an infinite molecular system, where the traction in the far-field can he intro¬ 
duced. This is motivated hy the ohservation that molecular simulations are typically 
conducted within part of the entire sample, due to the heavy computational cost. Math¬ 
ematically, the extra degrees of freedom in the surrounding region can he eliminated 
hy solving the finite difference equations associated with the molecular statics model. 
This gives rise to a BC, which is expressed as an extrapolation of the displacement to the 
atoms outside the boundary, along with a shift vector, which depends on the traction in 
the far field. We further demonstrate that the typical approach in which external forces 
are directly applied at the boundary might be incompatible wdth these BCs, and that 
they can lead to ill-posed problems. 

The present approach allows one to simulate a material system with local defects 
under traction load, which mimics a surround elastic medium. Another potential ap¬ 
plication is to the domain decomposition (DD) method for solving a large-scale molec¬ 
ular system, where the problem is divided into sub-problems, each of which is asso¬ 
ciate with a sub-domain. In particular, the Dirichlet-Neumann method and Neumann- 
Neumann method (e.g., see |[^) offer a coupling strategy without creating overlapping 
regions. Our method can be implemented within the DD framework to fascilitate par¬ 
allelization. 

The paper is organized as follows: We first consider a one-dimensional system to 
demonstrate how the BC can be derived. We further analyze the stability of the result¬ 
ing boundary value problem and the continuum limit. As an application, and a demon¬ 
stration of such BCs in high dimensions, we consider a dislocation dipole problem in 
section]^ We close the paper by a summary and some discussions. 

2. A ONE-DIMENSIONAL EXAMPLE 

To better illustrate the idea, let us first consider a one-dimensional semi-infinite chain 
of atoms with undeformed position x, = i, where i g Z+ u {0}. We will also use the un¬ 
deformed position to label the atoms. The deformed position is denoted by yt with 
displacement Ui. We assume that the atomistic potential has next-nearest-neighbor 
pairwise interactions. Namely, the total energy in the bulk can be written as, 

E(^(D+2-yd + f^(yi+i-yi))- 

I>0 

The extension to more general potentials and higher dimensions will be discussed later. 

Intuitively, there are at least two ways to impose a traction at the boundary. For in¬ 
stance, one may apply forces, denoted by Tq and Ti, to the first two atoms, located at 
yo and yi, respectively. Alternatively, one can introduce two additional atoms outside 
the boundary, which in the present case, have current positions y _2 and y_i. These ad¬ 
ditional atoms will be referred to as ghost atoms, since they play different roles as the 
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atoms in the interior. By specifying y _2 and y-i (or u-z and u-i), one also creates a 
traction at the boundary These two methods are illustrated in Fig.[^ 

To Ti 
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Figure 1. A one-dimensional chain with next nearest-neighhor inter¬ 
actions. Top: Tractions, in the form of point forces, are applied to the 
two atoms at the left boundary; Bottom: two ghost atoms are introduced 
outside the boundary. 


Remark 1. At this point, both methods seem plausible, and it is not immediately clear 
how they are related to each other. Also notable is that there are two parameters available 
to specific one traction condition. This will be clarified in the next section. 

We illustrate possible BCs with the following numerical experiment: We consider the 
one-dimensional chain model with harmonic interactions. The force constants are Ki - 
1 and Kz - -0.2 for the nearest neighbor and next nearest neighbor interactions. A 
traction force needs to be applied at the left boundary, while the atoms at and xjv+i 
are fixed. We choose N - 20. Three traction BCs are tested. In the first case, we apply 
a unit force on the first atom, and in the second case, we apply the same force to the 
second atom. In the third test, we split the force among the first two atoms and 
i). For such a simple setup, one would anticipate that the corresponding continuum 
model has a simple solution which is given by a uniform deformation gradient. These 
results are shown in Fig. In all these cases, the solutions develop a boundary layer, 
and none of them is fully consistent with the continuum solution. As comparison, we 
include the result from the traction BC that will be derived later, in which the position 
of the first two atoms is determined based on the given traction. It is clear that the 
boundary layer has been eliminated. 

2.1. The derivation of the traction BC. To understand the traction BC, we start by em¬ 
bedding the semi-infinite atom chain into an infinite chain. We recall that Xi - i, i eZ 
denotes the equilibrium positions. We take the view point that the BCs acting on the 
atom at Xq should be determined by the interaction of the semi-infinite chain with the 
atoms on the left (i < 0), whose degrees of freedom will be implicitly incorporated. In 
other words, the atoms Xi, i g Z_ serves as an environment for the system we consider. 
We will hence distinguish the two groups of atoms by referring them as system atoms 
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Figure 2. Testing several boundary conditions for the one-dimensional 
model. Left to right: A unit force applied to the first atom; A unit force 
applied to the second atoms; Forces applied to both atoms; Traction 
boundary condition that will be derived in the next section. 


and environment atoms, respectively. This is based on our observation that most atom¬ 
istic simulations are focused on part of the entire sample due to the small spatial scale 
associated with molecular models. 

The problem now is reduced to removing the atoms in the environment. To facilitate 
the reduction of degrees of freedom from the whole chain to the semi-infinite, we will 
take the harmonic approximation. This amounts to assuming that the interaction be¬ 
tween the atoms xi and Xj is harmonic if either / < 0 or j < 0. In this paper, we focus our 
attention to static problems, which can be formulated as an energy minimization prob¬ 
lem. For the current problem, the potential energy will be divided into several terms in 
accordance with the partition of the system, 

(1) E — iJsys + f^int + f^env 

where E^ys is the interaction among the atoms in the semi-infinite chain on the right: 

(2) Esys = iyi +2 - yi) +(yi+i - yd)- 

i>0 

Meanwhile Eint collects the interaction terms between a system atom and an environ¬ 
mental atom. In terms of the displacement uj - yj - xj, we have, 

2 ^2 2 ^2 2 
£’int= yiwo-M-l) +y(Wo-W- 2) -ty(Ui-U-i) , 


( 3 ) 
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where ki - y"(l) and jc 2 = V"[2) are two stiffness constants for nearest-neighbor and 
next-nearest-neighbor interactions. Further, ijenv denotes the energy for the environ¬ 
ment, given by 

(4) -Fenv ~ y^. ~ Uj-i) +—{Ui — Ui-2) ]• 

i<0 4 2 

Since the interaction is assumed to be next-nearest-neighbor, the environment atoms 
only interact with atoms with reference positions Xq and Xi, but not other system atoms. 
Given the displacement uq and ui, the force balance equations for the environment 
atoms can be written as 


(5) K2{.Uj+2-2Uj + Uj-2)+Ki{Uj+i-2Uj + Uj-i)-0, J G Z-. 

The general solution of this finite difference equation is given by, 

(6) Uj = A. + B j + CX^ + DX ^, 

where A is a root of the characteristic polynomial associated to ^ (the other roots are 
1 / A and 1 with multiplicity 2): 


(7) 




We collect some elementary properties of A in the following lemma. 


Lemma 1. Assume /ci > 0 and k -ki + 4x2 > have I A| < 1 and -K 2 A > 0. 


Proof. First notice that 


1-F 


4X2 , 4X2 4x^ I 2X2V 

-< l-F- + ^ =1-1--> 

Xi Xi xf '■ Xi 


which yields 



2X2 


For X 2 > 0, we then get 



Hence, by definition of A in 0 , we get A e [ -1,0]. 
For X 2 < 0, we have 
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This yields A > 0. Since 1 + 4 k2/ki < 1, we also have 


/ 4K2 4K2 

'l + — > 1 + — 


Kl 


T<1 


and hence 



We conclude that A < 1. 



□ 


Since | A| < 1, A^ grows exponentially as j -oo and hence unphysical. This leads to 
the requirement that C = 0. The positions of the atoms at Xq and Xi further provide two 
BCs for (0. The remaining one degree of freedom is determined hy the traction at the 
boundary 

(8) - Kiiuo - U-i) - Kziuo - U-z) - KziUi - U-i) = T, 

where T is the prescribed traction at the boundary (scalar in ID). Intuitively, the trac¬ 
tion across a material interface is given by the sum of the forces between two atoms 
that are on different sides of the interface |[^|^ . This is indeed non-trivial, especially 
for multi-body interactions. But formulas are available for most empirical potentials 

( 21 ). 


Remark 2. We further remark that the traction is conserved since no external force acts 
on the fictitious atoms: For j g Z_, 


^-KiiUj-Uj-i)-KziUj - Uj-z)^-KziUj+i - Uj-f 

^ - Uj) - KziUj+z- Uj)^ - KziUj+i - Uj-i) 

= -KiiUj+i - Uj)-Kz{Uj+i - Uj-i)-Kz{Uj+z-Uj). 

The solution to 10 can now be found. In particular, the coefficients are 


-A ^Uq + Ui + T/k 




B = - T Ik, and C = 


Uq-Ui-T/k 


1-A-i . , ^ ^ ^ , 

where k = ki+ 4kz. As a result, the displacements U-i and U-z are given by 
(9a) u_i = (l + A)uo-AMi + (l-A)r/K; 

(9b) M-2 = (1 + A)m-i - Amo + (1 - A)r/7C. 


Notice that U-i and U-z depend linearly on the displacements Uq, Ui and the traction 
T. This is a result of the harmonic approximation in the environment. Now that the 
degrees of freedom associated with the atoms further on the left are removed, we can 
formulate the boundary value problem for the semi-infinite atom chain xt, i eZ+u {0} 
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in terms of ghost atoms x-i and x -2 at the boundary. In general, the number of needed 
layers of atoms is determined by the interaction range. 

With the BCs, the molecular statics model is complete. We consider a slightly more 
general problem by allowing body forces fi to be applied to the system atoms. In this 
case, the force balance equations read as follows, 

(10a) - V’[yj+2 - yj) - V'iyj+i - yj) + V'iyj - yj-i) + V'iyj - yy_ 2 ) = fj, 7^2 

(10b) - V'{y 2 - yi) - l^'(y 2 - yi) + l^'(yi - yo) + - u-i) = /i, 

(10c) - l^'(y 2 -yo) - V'{yi-yo)+-Ki{uo-u-i)+K 2 {uo-u- 2 ) =/o, 

together with the BCs given by 

Another observation is that due to the semi-infinite nature, l[^-l[9| can only deter¬ 
mine u up to a constant. To uniquely fix the arbitrary constant, we choose Uq = 0. In 
addition, while the solution u can be a linear function that corresponds to a uniformly 
stretched (or compressed) state, it is natural to exclude those solutions that grow su- 
perlinearly at infinity |[^ . Hence, the complete set of BCs for the semi-infinite chain 
consists of the traction BC @ and the conditions 


(11a) 


uo = 0; 


(11b) 


limsup 

j^oo 



< OO. 


We emphasize that the above two BCs are associated with the semi-infinite chain 
under consideration, but not the traction at the left end of the chain. For finite system 
with a right boundary, appropriate BCs should be chosen to replace ([^ according to 
the physical situation. Our emphasis, however, is on the traction condition at the left 
boundary. 

Let us summarize the general framework for our BC construction in one-dimensional 
systems as follows 

Step 1. Supplement the system with a fictitious environment of atoms with linear ap¬ 
proximation; 

Step 2. Solve the positions of the environmental atoms with the condition of fixed trac¬ 
tion; 

Step 3. The BC of the atomistic system is then given in terms of the positions of the 
ghost atoms. 

This procedure can be clearly generalized to one-dimensional atomistic systems with 
arbitrary short-range interactions. The number of BCs depends on the interaction range. 

Next, we turn to several properties of the BCs. These will help us better understand 
the traction BC and also facilitate the extension to higher dimension. 
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2.2. The continuum limit. For continuum elasticity models, traction BCs are imposed 
as the normal component of the stress. In this subsection, we show that the continuum 
limit of the reduced system |[T^, together with the BCs l|9|, leads to the Cauchy-Born 
elasticity with continuum traction BC in elasticity. Hence, our BCs m can be viewed as 
the molecular statics analog of the traction boundary condition in continuum elasticity. 

To this end, we adopt the natural rescaling of the system such that the distance be¬ 
tween nearest-neighbor atoms in equilibrium becomes e. We will use superscript to 
make explicit the dependence on the scaling parameter e. Hence, the equilibrium po¬ 
sitions scale to = je, j eZ+u {0} and the deformed positions are y‘: = x‘j + u'^(xp. We 
rewrite the force balance equation and the traction BC accordingly: 



with 


(13a) {1 +X)Uq-X ul + e{l-X)T/ k; 

(13b) u^2 = (1 + - Auq + e(l - X)TIk. 

We now take the continuum limit £ —>■ 0. To the leading order, the equation ( |12al l be¬ 
comes 


(14) 


-div[y'(/(x)) + 2y'(2/(x))] =/(x). 


Note that for the current atomistic interaction potential, the Cauchy-Born stored energy 
density is given by : 

(15) Wcb(A] = V[! + A} + V[2(I + A]). 

Hence, ([T^ is exactly the force balance equation for the Cauchy-Born elasticity since 

(16) dAWcB(u') = ^'(1 + u'M) + 2V'(2 + 2u'(x)). 

Combining | |12bt , ( |12c| l and the BCs we get 


(17) 




rdyl-yl 




To the leading order, this yields 


(18) - ^'(1 + u'(0)) - 2V'[2 + 2u'(0)) = T. 

As the left hand side is equal to n ■ |^^o’ where n is the unit exterior normal, 

the BC is exactly the traction BC for the elastic energy density Wcb- 
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2.3. Linear stability at the equilibrium. We also make the observe that the force bal¬ 
ance equations l[^ with the traction BC can be viewed as a finite difference sys¬ 
tem with BCs. It is then natural to analyze the stability of such a finite difference sys¬ 
tem, similar in spirit to the analysis in the context of atomistic-to-continuum methods 
Hn). We also note that the stability is the crucial ingredient for the rigorous proof of 
the continuum limit in §2.2 The stability of molecular statics models under periodic 
and Dirichlet BCs have been analyzed in . 

To understand the stability, we linearize the force balance equations ([T^ at the equi¬ 
librium (undeformed) state, yielding, 

(19) -K2[Uj+2-2Uj + Uj-2)-Ki[Uj+i-2Uj + Uj-i)^ fj, j>0, 

supplemented by the BC @ and l[^ . Given T, we define the map Ht : (N) ^ (N) as 

(20) iHTU)j^-K2iUj+2-2Uj + Uj-2)-KiiUj+i-2Uj + Uj-i), j>0 

with u-i and u -2 determined by l[^ (and hence the dependence on T). Thus we have 

-(?C2(2-i-A)- l-/ci)(l - A)r/K, 7 = 0; 

-K2(l-A)r/?c, 7 = 1; 

0 , otherwise. 

Let us introduce a few short-hand notations. First we define the forward difference 
as 

(22) {Du)j-Uj+i-Uj. 

Moreover, the discrete Laplacian is given by 

(23) (^df ^)7 ~ u2uj Uj-i- 
Direct calculations yield, 

(24) {Ad2^dU )7 + 4(AdU )7 = ( 1 / 7 + 2 - 41 / 7+1 + 6 ^ 7 - 4 ^ 7-1 + 1 / 7 - 2 ) + 4 ( 1 / 7 +! - 21/7 + 1 / 7 - 1 ) 

= {Uj+2-2Uj + Uj- 2 ). 

With these preparations, we calculate the quadratic form {u,Hq{u))\ 

00 00 

{U,HqU) = -K2 Y, Uj[{^d^dU)j+^{^dU)j)-Kl Y Uj[Au)j 


( 21 ) 


Ht{u) j - Hoiu) I 


(25) 


1=0 


1=0 


= -K2 Y UjiAAu)j -k Y Uj(Au)j 
1=0 1=0 

Summation by parts gives (recall that Uq = 0) 

00 00 

^//7(A//)7 = -^|(D//)7|2. 
1=0 J=o 


( 26 ) 
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For the fourth order term, we get 

OO OO OO 

(27) Uj(AAu)j = Y - U-iiui + U-i) - ^ |(Au)y|^ + A(1 - X)ul, 

7=0 7=0 7=0 

where in the last equality, we used U-i = -Xui in the case T = 0 and Uq - 0. 

Finally, we have, 

OO OO 

(28) {u,Hou) =K^|(Du)j|^-K 2 ^|(Au)y|^-JC 2 A(l-A)Ui. 

7=0 7=0 

By Lemma[^ we have - 7 C 2 A(1 - A) > 0 as long as JCi > 0 and 7 c> 0. Thus, we have 

OO OO 

(29) {u,Hqu) >K^|(Du)j|^-7C2 ^|(Au)j|^ 

7=0 7=0 

Therefore, the scheme is stable as long as the underlying atomistic model is stable. 

For a general traction T, we have 

(30) {u,HTU)-{u,Hou)-—il-X)Tui. 

K 

Therefore, the stability follows from the stability in the case of T = 0. 

Remark 3. This analysis also shows that ifX in l[9| is replaced by an appropriate approx¬ 
imation, i.e., A a A satisfying-K 2 X{\ - A) > 0, a stable model would also be obtained. 


We show that a careless choice of the BC may lead to instability of the scheme. In¬ 
stead of let us consider an alternative set of BCs (to distinguish, we use ii for the 
displacement) 

(31a) U-i = (1 -F A “^)Uo - A “^Ui -t (A”^ -l)T /jc; 

(31b) U-2 = (1 + A“^) U-i - A“^ Uo -F (A”^ -l)T Ik. 

It is straightforward to check that this set of BC also yields traction T at the boundary 

(32) - Kiiuo - ii-i) - K 2 {uo - U- 2 ) - K 2 {ui - U-i) - T, 


and hence consistent with the traction BC in continuum elasticity. However, the result¬ 
ing scheme with the BC is not stable. In fact, we even lose uniqueness: it is easy to 
check that uj = A^ - 1 for j > -2 satisfies 


7 ^ 0 , 


(33) 


-K2iUj+2-2Uj -F Uj- 2 ) -KiiUj+i - 2Uj -F Uj-i) - 0, 
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and also at the boundary 

(34) S_i = (1 + A“^)uo - 

(35) u_2 = (1 + A”^)u_i - A”^Uo; 

(36) Uq - 0; 

\Uj\ 

(37) limsup —= 0. 

7—00 ] 

2.4. Connection to BCs with applied forces at the boundary. As we alluded to at the 
beginning of this section, it is also possible to apply forces (To and Ti) directly at the 
boundary to create a traction. But it is not immediately clear how much forces to apply 
on each of the two atoms at the boundary. Here we will demonstrate the connection to 
that approach. In particular, this discussion will also shed light on the selection of the 
forces. 

If we substitute l|9| into the force balance equation |[T^, we get 

-V'{y3 - yi) - I^'(y 2 - yi) + V'iy^ - yo) = /i + Ti; 

-v'iyz - yo) - I^'(yi - yo) = fo + n. 

with 

(38) Ti - K2iu-i - Uq) - K2((1 + A)(uo - Ui) + (1 - X)TIk]-, 

(39) To = K2{U-2-Uq)+Ki{U-i-Uq) 

= X{ki + K2(1 + A)) {uq — Ui) + (1 — A) + jc2(2 + A)) T/k. 

This provides the formulas for the forces. An important observation, however, is that 
these forces should depend on the displacement of the atoms at Xo and Xi. 

2.5. Traction BC from the Green’s function. To facilitate the extension of the BC to 
two dimensional systems, we take yet another point of view of the traction BC, from the 
lattice Green’s function perspective. 

Let us define the lattice Green’s function associated with the model 10, 

(40) -K2{Gj+2-2Gj + Gj- 2 )-KiiCj+i-ZCj + Gj-i) - 6{j), j 

In general, the lattice Green’s functions are useful analytical tools for studying lattice 
distortions around defects (see e.g., 1^). A typical route to compute the Green’s func¬ 
tion is via a Fourier transform, 

G(0 = E 

jeZ 


( 41 ) 
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with the inverse given hy, 
(42) 

This leads to 


Gj- 


I nn 

2,71 J-ji 




(43) 
and 

(44) 


- K2[e^^^ -2 + e“2'^)G(a -Ki- 2 + e-'^)G(a = 1, 


G(0 = 


4/C2sin^(0 +47Cisin^(^/2) 

However, due to the singularity at ^ = 0, the integral with G given above is not 
well defined. A remedy is to modify 


(45) 


G, = — 
^ 2n 


1 

■n j-7i 

1 

271 J-n 


-i(j 


47C2sin^(0 +47Cisin^(^/2) 
2 sin^ (0'/2) 


d^ 

d^ 


271 J-n 4jc 2 sin^(0 + 4 ki sin^ (^ 12) 

Conceptually, the Green’s function |[4^ is only defined up to a constant, and one can fix 
Go = 0 by subtracting a (infinite) constant from ( |42] l. 

As a result, the integral is now well defined as the integrand is regular as ^ ^ 0. The 
function Gj defined this way still satisfies the equations l [40t . We now make the con¬ 
nections to the BCs (|^. 

Lemma 2. For j < 0, 

(46) Gj-i = (1-t A)Gj - AGj+i 

Proof. Rewrite using a change of variable z = exp(i^) and the characteristic poly¬ 
nomial associated to the denominator, we get 

{z~-i - l)z 


Gj=—( — 

27TK21 J\z\ = l[Z- 

^lim f 

27lK2i e^OJr, 


1)2(z-A)(z-A-1) 
(z“-( - l)z 


dz 


■dz 


(z-l)2(z-A)(z-A-i) 

where the contour je is given by the boundary of Bi (0) \Be (1) on the complex plane. The 
second equality uses the fact that the integrand is regular as z ^ 1. 

Using this representation, we have 

z”''(z^-(l-t A)z-tA) 

-iiiii I 

271X2 i 


Gj-i — (1 -t A)Gj -t AG 


7+1 


— lim f 
TK2i 

— lim f 

K2i £^0jy 


(z-l)2(z-A)(z-A-i) 


dz. 


dz 


(z-l)(z-A-i) 
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As j <0 and |A| < 1, the integrand is holomorphic in Bi(0)\5e(l) for any e, and hence 
the integral vanishes for any e hy Cauchy’s integral theorem. Therefore, 1(4^ holds. □ 

The equation ( |46l l is exactly in the same form as the BCs @ when T -0. This is not 
surprising, since the Green’s function represents a special set of solutions. In particular, 
Gj satisfies the homogeneous difference equations Nevertheless, this simple ob¬ 
servation can he employed to determine the coefficients in the BCs hy using the Green’s 
functions as test functions. This will he implemented for problems in two-dimensions, 
and the implementation will be discussed in the next sections. 

3. Implementation in two-dimensional models 

Here we demonstrate how the BC can be extended to two-dimensional systems. 

3.1. The traction BC and the induced boundary value problem. For multi-dimensional 
problems, the BC is typically non-local |[^, in that the displacement of all the atoms at 
the boundary is coupled. It is also possible to consider nonlocal boundary conditions, 
for example, in the spirit of the boundary element method for molecular static models 
by one of the authors |j^ . Another alternative is to seek a local BC, in the sense that the 
position of the ghost atoms are only determined by the positions of nearby atoms in 
the system. To make the dependence local, we would employ a “local flattening’’ of the 
boundary. Roughly speaking, for an atom at the boundary, the position is determined 
by a homogeneous approximation of the local atom configuration with the local value 
of the traction tensor. 

To better explain the idea, we consider the face-centered cubic (FCC) lattice of Alu¬ 
minum with the axis aligned in (110>, <001> and <liO> orientations. When projected to 
the <liO> plane, the lattice spacing in the horizontal and vertical directions are and 

Y, respectively, which makes it look like a triangular lattice, as shown in Fig. Again, 
we introduce ghost atoms outside the boundary in order to achieve the desiredtraction 
condition. They are represented by open circles in Fig. 

Our main goal is to determine the actual position of the ghost atoms based on the 
displacement of the atoms in the interior and the applied traction T, which is a two- 
dimensional vector. In this case, it is in general cumbersome to obtain the exact bound¬ 
ary condition. Motivated by the one-dimensional traction BC, we seek an approximate 
BC in the following form, 

(47) Uj=Y,BjiUi + pj. 

ieSj 

The shift vector p is similar to the non-homogeneous term in l|9|, and it will be de¬ 
termined so that the correct traction is obtained. In the case when p-0, this boundary 
condition would coincide with the BCs that models an environment that is at a me¬ 
chanical equilibrium with zero stress 0I!3- In principle, an exact BC in this form can 
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Figure 3. The projected atomic position of a FCC lattice. Filled circles: 
Atoms in the interior; Open circles: the ghost atoms introduced outside 
the boundary. The boxes contain the set of atoms Sj that will be used to 
determine the displacement of the jth atom (see equation l[47t). 


be derived, e.g., in 114 , which mathematically, is a discrete analogue of the Dirichlet- 


to-Neumann (DtN) map. The exact expression is typically nonlocal, in that the summa¬ 
tion is over all the atoms near the boundary. But here we choose a local approximation, 
and restrict the summation in |[47) to those atoms that are close to the j atoms. These 
neighbors are collected by the set Sj. Due to the translational symmetry of the lattice, 
we will use the same set of neighbors when implementing the formula ( |47t , which is 
also demonstrated in Fig. More specifically, we start with the layer of ghost atoms 
closest to the boundary and apply the BCs ^7^. Once the displacement of these atoms 
are updated, we move to the next layer, and these steps will be repeated until the posi¬ 
tion of all the ghost atoms are updated. 

We now discuss how to determine the coefficients Bji. Since they are independent 
of the applied traction, they can be computed for the case T = 0 and p = 0. In this case, 
these coefficients can be determined using an optimization procedure, developed in 
0. More specifically, we choose an objective function as follows. 


(48) 


minfi, = 

k i 


Here Gj^k is the two-dimensional lattice Green’s function [18|. The main observation is 


that the BC should by satisfied by special solutions, especially the lattice Green’s func¬ 
tions Gj,fc, which corresponds to the solution of the linearized molecular statics model 
when a point force is applied on the A:th atom. This was already observed for the one¬ 
dimensional model. Ideally, the objective function would be zero when the BC is exact. 
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Further, we introduce a constraint, 

(49) 

i 

so that the constant solutions are admitted. This is also seen in the one-dimensional 
system: The two coefficients in @ add up to 1. 

It remains to estimate the vector p. In principle, it should he determined hy requiring 
the traction to arrive at a prescribed value. The total tractions along the boundary is 
given by the sum of the forces fT|[2T), 

(50) X fa- 


Here, fij comes from a force decomposition. Namely, 

(51) — — 

i^i 


This formula, which was already indicated by @, is consistent with the intuition of 
Cauchy. The explicit expressions for the force decomposition | [^ , especially for multi¬ 
body potentials, can be found in (2)[2T 


To control the local traction, we divide the region with ghost atoms into blocks, each 
denoted by a - 1,2,• • • ,M. The computational domain is denoted by n, and the 
intersection with is written as dn^. For each block we introduce a shift vector 
Pa- They are chosen so that the traction along dCa agrees with a prescribed value ta- 
This arrangement is illustrated in Fig. 


Ha-l 

1 

lift 

1 Ha+l 

1 

1 


1 


dQ(x 

3Qa+l 


Q 

Figure 4. Imposing tractions along the boundary: n indicate the com¬ 
putational domain. The atoms outside the boundary are grouped into 
blocks with the intersection with the boundary given by dDa- The 
traction on each dQa is prescribed. 
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We now put the mathematical models together. 


(52) 


d 

—y = o, 

dUi 

^ Uj— ^ji P a> ^ dUg;, 

i 

^ fij - ta- 

ieQyjeO,^ 


The first set of equations represent the force balance in the interior, with potential en¬ 
ergy given hy V. The remaining equations serve as BCs with prescribed tractions ta- The 
unknowns are the atomic displacement, together with the shift vectors pa- The atomic 
degrees of freedom associated with the atoms outside Q. has been implicitly taken into 
account by the second and third equations in In the next section, we will discuss 
an implementation method. 


3.2. Numerical implementations. Our reduced model | |52t consists of a set of nonlin¬ 
ear algebraic equations. It is therefore natural to make use of iterative methods, such as 
the quasi-Newton’s method. In general, this requires the approximation of the Jacobian 
matrix, since the analytical form is usually not available. The convergence is typically 
slow, especially when the system is not well prepared. 

To find an alternative, we notice that in the domain O, the molecular statics model 
is associated with an energy. Therefore, for Dirichlet BCs, where the atoms outside the 
boundary are held fixed, the solution to the molecular statics model correspond to an 
energy minimization, which is usually more robust and much more efficient than solv¬ 
ing the nonlinear equations. 

We implement the equations by a domain decomposition approach, and alternate 
among the three sets of equations in (j^. As an example to explain the idea, we may 
consider the coupling of the first two sets of equations and assume that pa is given. 
We create a few overlapping layers, in which the atoms serve as both ghost atoms and 
interior atoms. This is illustrated by Fig. In each iteration, we first update the dis¬ 
placement of all the ghost atoms including those in the overlapping region using 
(or the second equation in ij^). We then turn to the interior atoms, assuming that 
other atoms are held fixed (open circles in Fig. [^. By minimizing the energy, we obtain 
the updated position of the interior atoms. The numerical implementation has been 
done with the BFGS package j^. This iteration can be continued until convergence is 
reached. This is simply the Schwartz iteration (^) between the two models. 


3.3. Results from numerical experiments. As a test problem, we consider a disloca¬ 
tion dipole under a shear load. The atoms around the two dislocations with opposite 
Burgers vectors are shown in Fig. The embedded atom model (EAM) ||^ has been 
used as the interatomic potential. 
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Figure 5. An overlapping domain decomposition. 



Figure 6 . The atoms near the two dislocations. 

We first manually control the vector p-ipi,0), and observe the infiuence on the trac¬ 
tion changes. As a quasi-static loading, we increase pi with small increment and then 
solve the molecular statics model using the domain decomposition method described 
in the previous section. For each step, we also compute the traction at the boundary. 
The history of the total boundary traction is shown in Fig. We observe that the trac¬ 
tion increases as pi increases. Flowever, when pi reaches certain critical value, a sud¬ 
den drop is observed. In this case, the two dislocations move to the boundary, and the 
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entire sample undergoes a complete slip. Fig. [^shows the atomic positions before and 
after the slip. 



Figure 7. The history of the traction at the boundary. 



Figure 8 . The position of the atoms before and after the yield stress. 

In the next experiment, we apply a uniform traction along the upper and lower bound¬ 
aries. In Fig. 1^ we show the resulting values of pi along the boundary. While, the result¬ 
ing tractions have reached the prescribed values (2 x 10 “^, 0), it is clear that the values 
for p are not homogeneous, mainly due to the presence of the dislocations. The dis¬ 
placement is shown in Fig. together with a close-up view of the atomic positions. 
All these results suggest that the atomic positions are not uniform. Compared to the 
simulations of dislocation dipoles using the Parrinello-Rahman method (e.g.,|j^), the 
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current approach does not introduce periodic images of the dislocation dipole. More¬ 
over, the uniform shear stress can he applied without forcing a uniform deformation 
along the boundaries. 



X 




Figure 9. The traction (top) and the shift vector pi (bottom) along the 
upper boundary. 


4. Summary AND Discussion 

We have formulated boundary conditions that impose a traction force on a molecular 
statics model. These boundary conditions are derived by taking into account the sur¬ 
rounding elastic environment. Flence, the computational domain is part of a much big¬ 
ger sample, and artificial boundary effects can be eliminated. In the continuum limit, 
these boundary conditions coincide with the Neumann boundary condition for contin¬ 
uum elasticity models. We have restricted our discussions to static problems. Extension 
to dynamic problems will be considered in future works. 
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Figure 10. The displacement ui and a close-up view of the atoms near 
the two dislocations (generated in Ovito 
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